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ABSTRACT 

We present a new approach for image reconstruction and weak lensing measurements with interfer- 
ometers. Based on the shapelet formalism presented in Refregier (2001), object images are decomposed 
into orthonormal Hermite basis functions. The shapelet coefficients of a collection of sources are si- 
multaneously fit on the uv plane, the Fourier transform of the sky brightness distribution observed by 
interferometers. The resulting x^"fit is linear in its parameters and can thus be performed efficiently by 
simple matrix multiplications. We show how the complex effects of bandwidth smearing, time averaging 
and non-coplanarity of the array can be easily and fully corrected for in our method. Optimal image 
reconstruction, co-addition, astrometry, and photometry can all be achieved using weighted sums of the 
derived coefficients. As an example we consider the observing conditions of the FIRST radio survey 
(Becker, White & Helfand 1995; White et al. 1997). We find that our method accurately recovers the 
shapes of simulated images even for the sparse uv sampling of this snapshot survey. Using one of the 
FIRST pointings, we find our method compares well with CLEAN, the commonly used method for in- 
terferometric imaging. Our method has the advantage of being linear in the fit parameters, of fitting all 
sources simultaneously, and of providing the full covariance matrix of the coefficients, which allows us 
to quantify the errors and cross-talk in image shapes. It is therefore well-suited for quantitative shape 
measurements which require high-precision. In particular, we show how our method can be combined 
with the results of Refregier & Bacon (2001) to provide an accurate measurement of weak lensing from 
interferometric data. 

Subject headings: cosmology: large-scale structure of universe - gravitational lensing - methods: data 
analysis - techniques: interferometric 



1. INTRODUCTION 

Interferometers are widely used for astronomical ob- 
servations as they provide high angular resolutions and 
large collecting areas. Existing interferometers in the radio 
(e.g. the Very Large Array (VLA), the Berkeley-IUinois- 
Maryland Association (BIMA), etc.) and in the optical 
band (e.g., the Center for High Angular Resolution As- 
tronomy (CHARA), the Palomar Testbed Interferometer, 
the Keck Interferometer, the Very Large Telescope (VLT) 
Interferometer, etc.) will soon be complemented by new 
facilities such as the Expanded Very Large Array (EVLA) , 
the Square Kilometer Array (SKA), the Atacama Large 
Millimeter Array (ALMA) , and the Low Frequency Array 
(LOFAR). Interferometers are now also being developed 
to produce maps of the Cosmic Microwave Background 
on small scales (e.g., the Cosmic Background Imager 
(CBI), the Array for Microwave Background Anisotropy 
(AMiBA), the Very Small Array (VSA), etc). 

Interferometric arrays, however, do not provide a di- 
rect image of the observed sky, but instead measure its 
Fourier transform at a finite number of discrete samplings, 
or 'uu' points, corresponding to each antenna pair in the 
array. The image in real space must therefore be recon- 
structed from the uv plane by inverse Fourier transform 
while deconvolving the effective beam arising from the fi- 
nite sampling (see Thompson et al. 1986, Perley et al. 
1989, and Taylor et al. 1999 for reviews). For this pur- 
pose several elaborate methods have been developed. For 



instance, the commonly used CLEANing algorithm imple- 
mented in the NRAO AIPS software package (Hogbom 
1974; Schwarz 1978; Clark 1980; Cornwell 1983), relies on 
successive subtraction of real-space delta functions from 
the uv plane. Another method is based on Maximum En- 
tropy (e.g., Cornwell & Evans 1985) and consists of find- 
ing the simplest image consistent with the uv data. These 
methods are well-tested and appropriate for various ap- 
plications; however, the methods are non-linear and do 
not necessarily converge in a well-defined manner. Con- 
sequently, they are not well-suited for quantitative image 
shape measurements requiring high precision. In particu- 
lar, weak gravitational lensing (see Mellier 1999; Bartel- 
mann & Schneider 2000 for reviews) requires the statistical 
measurements of weak distortions in the shapes of back- 
ground objects and thus cannot afford the instabilities and 
potential biases inherent in these methods. While inter- 
ferometric surveys offer great promises for weak lensing 
(Kamionkowski et al. 1998; Refregier et al 1998; Schnei- 
der 1999), a different approach for shape measurements is 
therefore required to achieve the necessary accuracy and 
control of systematics. 

In this paper, we present a new method for reconstruct- 
ing images from interferometric observations. It is based 
on the formalism introduced by Refregier (2001, Paper I) 
and Refregier & Bacon (2001, Paper II), in which object 
shapes are decomposed into orthonormal shape compo- 
nents, or 'shapelets'. The Hermite basis functions used 



in this approach have a number of remarkable properties 
which greatly facilitate the modeling of object shapes. In 
particular, they are invariant under Fourier transformation 
(up to a rescaling) and are thus a natural choice for inter- 
ferometric imaging. We show how shapelet components 
can be directly fitted on the uv plane to reconstruct an in- 
terferometric image. The fit is linear in the shapelet coeffi- 
cients and can thus be performed by simple matrix multi- 
plications. Since the shapelet components of all sources are 
fitted simultaneously, cross-talk between different sources 
(e.g., when the sidelobe from one source falls at the posi- 
tion of a second source) arc avoided, or at least quantified. 
The method also provides the full covariance matrix of 
the shapelet coefficients, and is robust. We also show how 
the complex effects of bandwidth smearing, time averaging 
and non-coplanarity of the array can be easily and fully 
corrected for in our method. Our method is thus well- 
suited for applications requiring unbiased, high-precision 
measurements of object shapes. In particular, we show 
how the method can be combined with the results of Paper 
II to provide a clean measurement of weak gravitational 
lensing with interferometers. We test our methods using 
both observations from the FIRST radio survey (Becker 
et al. 1995; White et al. 1997) and numerical simulations 
corresponding to the observing conditions of that survey. 
We also show how our method can be implemented on par- 
allel computers and discuss its performance in comparison 
with the CLEANing method. 

Our paper is organized as follows. In §2, we first sum- 
marize the relevant features of the shapelet method. In §3, 
we describe how shapelets can be applied to image recon- 
struction with interferometers. In §4, we discuss tests of 
the method using both simulated and real FIRST observa- 
tions. In §5 we show how our method can be used for weak 
lensing applications. Our conclusions are summarized in 
§6. 

2. SHAPELET METHOD 

We begin by summarizing the relevant components of 
the shapelet method described in Paper I. In this ap- 
proach, the surface brightness /(x) of an object is decom- 
posed as 



/(x) = > /„i?„(x;/3). 



where 



5n(x;/3) 



7I„,(/3-ixi)i/„,(/3-ix2)e ^ 



(1) 



(2) 



are the two-dimensional orthonormal Hermite basis func- 
tions of characteristic scale /3, i?m(C) is the Hermite poly- 
nomial of order m, x = (xi, 0:2) and n = (ni, 712). The ba- 
sis is complete and yields fast convergence in the expansion 
if /3 and x = are, respectively, close to the size and loca- 
tion of the object. The basis functions can be thought of 
as perturbations around a two-dimensional Gaussian, and 
are thus natural bases for describing the shapes of most as- 
tronomical objects. They are also the cigcnfunctions of the 
Quantum Harmonic Oscillator (QHO), allowing us to use 
the powerful formalism developed for that problem. A sim- 
ilar decomposition scheme using Laguerre basis functions 
has been independently proposed by Bernstein & Jarvis 
(2001). 



The Hermite basis functions have remarkable math- 
ematical properties. In particular, let us consider 
the Fourier transform of an object intensity, /(k) = 
(27r)~2 J_ (Pxf{x)e^^'^. It can be decomposed as 

/(k) — X^n/n ^n(k;/3), whcrc i?„(k;/3) are the Fourier- 
transforms of the basis functions, which obey the dual 
property 

5„(k;/3)=z("i+"^)i?„(k;/3-i). (3) 

From the orthonomality of the basis functions, the coeffi- 
cients are given by 



/n 



d2fc/(k)B„(k;/3). 



(4) 



This invariance (up to a rescaling) under Fourier transfor- 
mation (Eq. [3]) makes this basis set a natural choice for 
interfcrometric imaging. 

3. SHAPELET RECONSTRUCTION WITH 
INTERFEROMETERS 

In this section, we describe how shapelets can be ap- 
plied to interfcrometric imaging. We first briefly discuss 
how images are mapped onto the uv plane by interfer- 
ometers. We also show how the uv plane can be binned 
into cells to reduce computation time and memory require- 
ments. We then describe how the shapelet coefficients can 
be directly fit onto the binned uv plane using a linear x^ 
procedure. Finally, we describe how the resulting shapelet 
coefficients can be optimally combined to reconstruct the 
image, to co-add several pointings, and to measure shape 
parameters. 

3.1. Interferometric Observations 

An interferometer consists of an array of antennae whose 
output signals are correlated to measure a complex 'visibil- 
ity' for each antenna pair (see Thompson et al. 1986, Per- 
ley et al. 1989, and Taylor et al. 1999 for reviews). Each 
visibility is then assigned a point on the 'uu plane' cor- 
responding to the two-dimensional spacings between the 
antennae. In practice, the visibilities are close to, but 
not exactly equal to a two-dimensional Fourier transform 
of the sky brightness. Within the conventions of Perley, 
Schwab & Bridle (1989) for the VLA, the visibility mea- 
sured for the antenna pair {i,j) at time t and at frequency 
ly is indeed given by 

(5) 
whcrc /(I, V, t) is the surface brightness of the sky at 
position 1 = (/, m) with respect to the phase center, 
and A(l, v) is the (frequency-dependent) primary beam. 
For the VLA, the primary beam power pattern can be 
well-approximated as the Bessel function 2Ji(z)/z, where 
z ~ 3.234 rd-^, 0p = 30'.83 x (^^^^), i^ is the observa- 
tion frequency and r is the position offset from the phase 
center (Condon et al. 1998). The u,v,w coordinates are 
given by 

(u\ I sin iJo cos Hq \ 

w = — sin (5o cos _ffo sin 8q sin Hq cos (5o 
w ) \ cos (5o cos Hq — cos b^ sin Hq sin (5o / 



(6) 




where A = cv^^ is the wavelength of observation, Hq and 
^0 are the hour angle and declination of the phase center, 
and {Lx, Ly,Lz) are the coordinate differences for the two 
antennas. The latter are measured in a fixed-Earth co- 
ordinate system, for which the sky rotates about the L^ 
axis. Note that the (u,w,w) positions of the visibilities 
define the synthesized beam pattern. Since the {u,v,w) 
coordinates are entirely determined by the antenna posi- 
tions, source coordinates, and time and frequency of the 
observations, the synthesized beam is precisely known for 
interforemeters. 

Only in the absence of a primary beam [A = 1), for ob- 
servations at zenith (w = 0), and for small displacements 
from the phase center {l,m <C 1), does the visibility reduce 
to an exact Fourier Transform of the intensity. Further- 
more, the visibilities are measured in practice by averaging 
over small time and frequency intervals. The resulting av- 
eraged visibility is given by 



Am) as 



V^,^ 



dt / dv T{t)G{v)V^j{t,v) 



(7) 



where T{t) and G{t) are the time and frequency window 
functions, respectively, and arc normalized as J dtT{t) = 
J di'G{i') = 1. Because the time and frequency intervals 
are typically very small, this double integral can be evalu- 
ated by Taylor expanding Vij {t, v) about the central values 
to and fg of the window functions. For square-hat window 
functions of width Ai (exact) and Ar^ (approximate), re- 
spectively, we obtain 



r/^2 



Vij ~ Vij{tQ,vo) + 



24 



52%(to,f^o) 



dt^ 



{Atf 






(Aiy) 



(8) 



When the telescope points to a fixed location on the sky, 
the hour angle of the phase center changes as Ho{t) ex Lust, 
where loe is the angular frequency of the Earth. On the 
other hand, the declination Sq of the phase center remains 

constant. 

The above expression for Vij can thus be computed an- 
alytically, leaving the two-dimensional 1-integral to eval- 
uate numerically. Note that this provides a direct and 
complete treatment of primary beam attenuation, time- 
averaging, bandwidth smearing and non-coplanarity of the 
array. These effects are difficult to correct for in the con- 
text of the standard CLEANing method. 

3.2. Binning in the uv plane 

In practice, the number of visibilities per observation is 
large (> 10^). Directly fitting the shape parameters to all 
uv points would thus require prohibitively large comput- 
ing time and memory. Instead, we use a binning scheme 
to reduce the effective number of uv points without loos- 
ing information. In the uv plane, we set a grid of size 
Am = i Al~^ and average the visibilities inside each cell, 
where A^ is one-half of the intended field of view, and the 
factor 2 accounts for the Nyquist frequency. The choice of 
Am is designed both to minimize the number of cells and to 
avoid smearing at large angular scales, which would oth- 
erwise act like an effective primary beam attenuation. We 
thus calculate the average visibility in the uv cell c (of size 



^-^Sn 



(9) 



ijec 



where Nc is the number of visibilities in the cell. This is 
the data we will use to reconstruct the image. 

3.3. Fitting for the shapelet coefficients 

We now wish to model the intensity fsilTin) of each 
source s as a sum of shapelet basis functions 

/s(l)=^./nsBn(l-ls;A), (10) 

n 

centered on the source centroid Is = (Is, ms), and scale (3s- 
Our goal is to estimate the shapelet coefficients /„« of the 
sources given the binned uv data {Vc}- (We will describe 
how the centroid and shapelet scales are chosen in practice 
in §4.1). In principle, the full mm plane provides complete 
shape information for the sources. However, due to the 
finite number and non-uniform spacings of the antennae, 
the uv (Fourier) space is poorly sampled, thus hampering 
the decomposition. This prevents us from performing a 
simple linear decomposition as is done with optical images 
in real space (see Paper I). This problem can be largely 
resolved by making a linear fit to the mm plane with the 
shapelet coefficients as the free parameters. 

For this purpose, the first step is to compute the binned 
visibilities V^ corresponding to each shapelet basis func- 
tions i?ns(I — ls;/3s) for each source s. This can be done by 
first computing the time- and frequency-averaged visibility 
V^j by setting /(I) — i?ns(l — Is; /3s) in Equations (5) and 
(7). To prevent potential biases introduced by the binning 
scheme, we evaluate the basis functions at every visibility 
point and then average them inside each cell to compute 
V^ just as in Equation (9). Note that this ensures that 
the systematic distortions induced by the primary beam, 
bandwidth smearing, time-averaging and non-coplanarity 
can all be fully corrected in our method. 

The next step is to form and minimize 

X^ = (d - M ff C-i (d - M_/), (11) 

where d = {Vc} is the data vector, M = {V^ } is the the- 
ory matrix, and f = {/ns} is the parameter vector. The 
covariance error matrix 

C = cov[d,d]=((d-(d))^(d~(d))) (12) 

for the binned visibilities can be estimated in practice ei- 
ther from the distribution of the visibilities in each bin 
or from the error tables provided by the interferometric 
hardware. 

Because the model is linear in the fitting parameters, 
the best-fit parameters f can be computed analytically as 
(e.g., Lupton 1993) 

f == (M^C-iM)-^M'^C-M. (13) 

The covariance error matrix W = cov[f , f] of the best-fit 
parameters is given by 

W'= (M'^C^^M)"^ (14) 

This provides us with an estimate for the shapelet coeffi- 
cients for each source and for their full covariance matrix. 
Note that all sources are fitted simultaneously thus avoid- 
ing (or at least quantifying) potential cross-talk between 
different sources (e.g., when a sidelobe from one source 
falls at the position of a second source). The coefficient 
covariance matrix can also be used to determine degenera- 
cies produce by the finite uv sampling of the array. 



3.4. Combining the Shapelet Coefficients 

Now that we have derived estimates /„ for the shapelet 
coefficients /„ for each source in a pointing, we can com- 
bine them to construct an image and to compute useful 
quantities such as the fluxes, centroids and sizes of the 
sources. 

We first consider the practical problem of co-adding sev- 
eral pointings to derive an optimal image of a source. Let 
fnp be the coefficients of a source derived from pointing 
p, and let W^mp be the associated covariance error matrix 
(from Eq. (14)). It is easy to see that the error in the co- 
added coefficients /„ will be minimized if they are given 
by the weighted sum 



./n = 






(15) 



nnp 



The covariance error matrix Wnm — cov[/„, /m] of the 
co-added coefficients are then given by 



W^n 



L^p * * nnp * ^ nmp ' » nimp 



(16) 



{EpWnrlp) (EpWn^L 

We can then find an optimal weighting to reconstruct 
the image of a source from the estimated coefficients /„. 
To do so we seek the reconstructed coefficients given by 

/n = <Pnfn. (17) 

The weights 0n are chosen so that the reconstructed im- 
age /''(I) ~ X^n /n^n(l) ^^ '^^ closc as possiblc' to the true 
image /(I), in the sense that the least-square difference 

fd'l [r(l)-/(l)f-E[/n-/n]' (18) 

•^ n 

is minimized. It is easy to show that this will be the case 
when 



l/nl 



l/nl 



Wr, 



l./n| 



Wn 



l/n 



(19) 



where the right-hand side provides an approximation 
which can be directly derived from the data. This weight- 
ing amounts to Wiener filtering in Shapelet space, in anal- 
ogy with that performed in Fourier Space (see, e.g.. Press 
et al. 1987). Figures 1 and 3 show several reconstructed 
images using this weighting scheme. Note that this pro- 
duces an estimate for the deconvolved image of the source. 
For display purposes, it is sometimes useful to smooth 
the reconstructed image by a Gaussian kernel (the 'restor- 
ing beam' in radio parlance). This can easily be done in 
shapelet space by multiplying the coefficients by the ana- 
lytic smoothing matrix described in Paper I. 

While Wiener filtering yields an optimal image recon- 
struction, it is not to be used to measure source parameters 
such as flux, centroid, size, etc. Instead, an unbiased esti- 
mator for shape parameters can be derived directly from 
the shapelet coefficients (see Paper I) . For instance, an es- 
timate for the flux F = J d^lf{V) of a source is given by 

F — A^f where 



A -7r5/325(2-»i-«2) f "■! 



«2 

n2/2 



, (20) 



if ni and n2 are both even (and vanishes otherwise). The 
variance uncertainty in the flux is then simply 

a^[F] ^A^WA, (21) 



which provides a robust estimate of the signal-to-noise 
SNR= F/a[F] of the source. Similar expressions can be 
used to compute the centroid and rms size of the source. 
This can be easily generalized to compute in addition the 
major and minor axes of the source and its position an- 
gle. Note that these expressions are, again, estimates for 
deconvolved quantities. 



4. TEST OF THE METHOD 

As an application, we consider the FIRST radio survey 
(Becker et al. 1995; White et al. 1997), being conducted 
with the VLA at 1.4 GHz in the B conflguration. For 
this survey, the primary beam FWHM is ^ 30' and the 
angular resolution is 5". 4 (FWHM). The survey currently 
contains about 7.2 x 10^ sources with a 5a flux-density 
limit of 1.0 mJy over A ~ 8,000 deg^; the mean source 
redshift is (z) ^ 1. Observing time has been allocated 
to extend its coverage to 9,000 deg^. The survey is com- 
posed of 165-second 'grid-pointings' with a time- averaging 
interval At = 5 seconds. It was conducted in the spec- 
tral synthesis mode, with a channel bandwidth of Av = 3 
MHz. Because this wide-field survey was performed in the 
snapshot mode, its uv sampling is very sparse. This makes 
shape reconstruction particularly challenging for FIRST, 
providing a good test for our method. 

As explained in §3.1, higher order effects such as band- 
width smearing and time-averaging produce small distor- 
tions in the reconstructed image shapes if they are left un- 
counted for. These must be carefully corrected for high- 
precision statistical measurements of object shapes such 
as those required in weak lensing surveys. The effects are, 
however, very small and not noticeable on an object-by- 
object basis. For the purpose of this test, we thus ignore 
these effects and instead focus on the dominant factor in 
shape reconstruction, the finite and discrete uv sampling. 



4.1. Simulations 

As a first test, we generated simulated VLA data using 
the observational parameters of FIRST. A grid pointing 
was generated at zenith with 33 5-second time-averaging 
intervals and 14 3-MHz channels in the B configuration. 
Simulated sources were randomly distributed within 23'. 5 
of the phase center, the cutoff adopted for creating the final 
co-added FIRST maps (Becker et al. 1995); the number 
density, flux density and size distributions chosen for the 
sources were similar to sources in the FIRST catalog. After 
generating the visibilities, we added uncorrelated Gaussian 
noise to the real and imaginary component of each uv data 
point, with an rms of a^ = (j„7V^jj?, where TVvis is the total 
number of visibilities. The real-space rms noise cr„ was set 
to 0.3 mJy beam^^, which is somewhat higher than the 
typical FIRST noise level, ~ 0.2 mJy bearn"^. 

We then simultaneously fitted all 23 sources in the 
grid pointing directly in the uv plane. We imposed the 
constraint that that source intensities are real (i.e., non- 
imaginary). Each source s was modeled as a shapelet with 
scale /3s, maximum shapelet order ninax,s7 and center po- 
sition Is,. In principle, it is possible to determine these 
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Fig. 1. — Three sources from the simulation of a FIRST grid pointing. The input images (before the addition of noise), dirty images, and 
shapelet reconstructions are shown from top to bottom, respectively. The images are 32" across with a 1" pixel size, and the resolution is ~ 
5". 4 (FWHM). The input source flux densities are (16.0, 4.3, 3.7) mJy from left to right, and the recovered shapelet flux densities calculated 
using Equation (20) are (16.1, 4.1, 3.8) mJy, respectively, where (10, 21, 6) shapelet coefficients were used in the reconstruction. The noise 
level in the simulation is about 0.3 mJy beam"^. 



parameters with a source detection algorithm which di- 
rectly uses shapelets. One can, for instance, tile ground- 
state shapelets with different sizes in the uv plane, and 
thus detect sources with different sizes. However, this is 
computationally expensive and, since the FIRST catalog 
is conveniently available, we have not implemented this 
algorithm. 

Instead, good choices for these parameters were derived 
from the FIRST catalog, which lists basic shape parame- 
ters for each source, such as its ccntroid, flux density, ma- 
jor and minor axes, and position angle, all obtained from 
an elliptical Gaussian fit. The shapelet position l^. was 
simply set to the centroid position from the catalog. The 
choices for the shapelet scales f3s and maximum shapelet 
orders nmax,s were derived as follows. As described in Pa- 
per I, the Hermite basis functions have two natural scales: 
^max corresponding to the overall extent of the basis func- 



tions, and ^niin corresponding to the smallest-scale oscilla- 
tions in the basis functions. These scales are related to the 
shapelet scale and maximum order by 6'niax ~ /3(^max + 1) 
and ^inin ^ /3(?^max + 1)^^- As rimax iucreascs, the large- 
scale size of the shapelet grows, while its small-scale fea- 
tures become finer. The shapelet thus becomes more ex- 
tended both in real and in Fourier space. We therefore 
choose ^max to be the rms major axis from the FIRST 
catalog, and ^min to correspond to the longest basehne of 
the VLA: ~ 1".8 (rms) in real space. This provides us with 
a choice for /3 ~ (6'inax Omin)°-^ and for rimax ^ f^ - 1 
for each source. 

Solving Equation (11), we then obtain the shapelet co- 
efficients and the covariance matrix using Equations (13) 
and (14). The results are presented in Fig. 1, where 
the input images (before the addition of noise), in- 
verse Fourier-transformed uv data ('dirty' images), and 



shapelet-reconstructed images (with Weiner filtering, see 
§3.4) of three of the sources are shown. Each image shown 
is 32" across and the resolution is about 5". 4 (FWHM). 
The poor uv sampling of FIRST and the effect of noise 
are evident in the dirty images. For both resolved (left 
panels) and unresolved or marginally resolved (right pan- 
els) sources, the reconstructions agree with the inputs very 
well. The more complicated structure in the central panel 
is not fully recovered by shapelets. This is expected, since 
the small-scale structure of the source is not resolved and 
therefore can not be fully recovered in the reconstruction. 
The comparison between the input and shapelet- 
reconstructed flux density for all sources in the grid point- 
ing is shown in Figure 2. The shapelet flux density is given 
by Equation (20) and its la error by Equation (21). The 
source flux densities are well-recovered by the shapelets in 
an unbiased manner. Note the range of error bars at a 
given input flux is due to the range of source sizes. For 
instance, for an input flux density of about 2 mJy, the 
source with a relatively large error bar has a major axis of 
about 8" (FWHM), while those with small error bars are 
unresolved or barely resolved (major axis FWHM < 5"). 
In general, we find the shapelet reconstruction from the 
sparsely sampled and noisy simulated data to be in good 
agreement with the input (noise-free) image. Note that 
our method can be used to identify and discard spurious 
sources arising from sidelobes and other artifacts in the 
dirty image. Indeed, when we place an extra shapelet cen- 
tered at a random positions in the field, the coefficients of 
that shapelet are consistent with zero. 



Flux Comparison 



o 



a; 
K 






Input, Flux (mJy) 

Fig. 2. — Flux density comparison for the simulations: the in- 
put flux densities of the 23 sources in the simulated pointing are 
compared with the shapelet reconstructed flux densities. The solid 
line corresponds to a perfect reconstruction- i.e., to a recovered flux 
density equal to the input flux density. 

4.2. Data 




including the primary beam response) of 0.75 mJy^. For 
each of the resulting 23 sources, we use the source major 
axis to estimate j3 and rimax as described in the previous 
section. We then simultaneously fit all the sources for the 
shapelet coefficients directly in the uv plane. Note that the 
shapelet coefficients obtained are deconvolved coefficients. 

Figure 3 shows the reconstruction of three representative 
sources in the bottom panels. Also shown for comparison 
are the images of the sources constructed using the stan- 
dard MPS CLEAN algorithm with a CLEANing limit of 
0.5 mJy (central panel), along with the dirty images (top 
panel). Each panel is 32" across and the FWHM of the 
FIRST resolution is 5". 4. The shapelet method does not 
involve image pixels in the modeling; one is therefore free 
to specify the pixel size when constructing the images. 
Here the dirty and CLEANed images have pixel sizes of 
1".8, while the shapelet images have pixel sizes of 1" and 
thus show finer details. For demonstration, the shapelet 
reconstructions have been Weiner-filtered using the result- 
ing covariance matrix. For a direct comparison, they have 
also been smoothed with a Gaussian kernel with a stan- 
dard deviation of 2". 3, reproducing the 'restoring beam' 
of the CLEANed image. We find that the shapelet recon- 
structions compare well with the CLEANed images. 

In further tests, we have encountered cases where a 
bright source (> 100 mJy) lies in or near a grid point- 
ing. We have found that the presence of the bright source 
does not affect the fit of the other sources in the grid in a 
noticeable way. Our method can thus well handle the dy- 
namical range of the FIRST survey, which spans more than 
3-orders of magnitude. For fainter sources (< 1 mJy; i.e., 
detection SNR < 6), the reconstructions are rather poor 
at times, in contrast to those for brighter sources (which 
are almost always well fitted). This is of course reasonable, 
given the larger impact of noise for faint sources. 

In Figure 4 we display a portion of the covariance ma- 
trix for the shapelet coefficients for the nine sources in the 
pointing with the highest peak flux densities. The hori- 
zontal and vertical lines separate the nine sources. The 
diagonal line from the lower-left to the upper-right corner 
represents the variance of the shapelet coefficients. The 
block-diagonal boxes are the covariance matrix of the co- 
efficients of the nine sources. The off-diagonal blocks quan- 
tify the cross-talk between sources. Note that the correla- 
tion between coefficients are roughly an order of magnitude 
smaller than the variance. Figure 5 shows the error in the 
shapelet coefficients (nl,n2) of the source shown in the left 
panels of Figure 3. (These errors are the diagonal segment 
of the 4*'* diagonal box counting from the lower left in 
Fig. 4). In general, we find that higher-n coefficients tend 
to be noisier. This is expected since convolution (or, equiv- 
alently, uv sampling) suppresses the small scale informa- 
tion encoded by coefficients with large n (see paper I) . The 
covariance matrix thus provides us with useful information 
on the error in each coefficient, and quantifies cross-talk 
between coefficients both within and among sources. 



Next, we test our method by applying it to one of the „ 

FIRST grid pointings (14195+38531). For this purpose, '^••^- '-computation 

we selected all sources within 23'. 5 of the phase center from Since the shapelet coefficients of all sources are simulta- 

the FIRST catalog with a measured flux density limit (i.e., neously flt to a large number of visibilities, the computing 

^ As explained in Becker et al. (1995), a map flux density of 0.75 mJy corresponds to a source flux density of 1.0 mJy owning to "CLEAN 
bias" corrections. 
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Fig. 3. — Three sources from one of the FIRST grid pointings. The dirty images, CLEAN images, and shapelet reconstructions are shown 
from top to bottom, respectively. The images are 32" across and the resolution is about 5". 4 (FWHM). The dirty and CLEAN images are 
displayed with a 1".8 pixel size, while the shapelet reconstruction images have 1" pixels. The source flux densities measured by integrating 
fitted elliptical Gaussians to CLEANed sources are (13.4, 14.2, 19.5) mJy, from left to right, and the recovered shapelet flux densities calculated 
using Eq. (20) are (17.0, 13.8, 22.5) mJy, respectively. We used (15, 28, 6) shapelet coefficients in the reconstructions, along with Wiener 
filtering and smoothing by a Gaussian restoring beam with a standard deviation of 2.3". 



memory required for the calculation is not negligible. We 
have implemented the method on the UK COSMOS SGI 
Origin 2000 supercomputer, which has 64 RIOOOO MIPS 
processors with a shared-memory structure. Numerically, 
the shapelet coefficients can be obtained by performing 
simple matrix operations as in Equation (13), or by solving 
the linear least-squares problem, M / = d, using matrix 
factorization or singular value decomposition, and assum- 
ing that the data covariance matrix C is diagonal. Both 
methods can be efficiently parallelized. With our binning 
scheme, the run-time memory required for this particular 
FIRST grid pointing was about 700 MB, for 23 sources 
and a total of 177 shapelet parameters. The CPU time re- 
quired was about 26 seconds with 10 processors or about 
5 minutes in scalar mode. For other grid pointings with 
different numbers of sources, the computation time ranges 



between 20 and 60 seconds with 10 processors, with a run 
time memory between 0.5 to 1.5 GB. 



5. APPLICATIONS TO WEAK LENSING 

Weak gravitational lensing is now established as a pow- 
erful method for mapping the distribution of the total mass 
in the Universe (for reviews see Mellier 1999; Bartelmann 
& Schneider 2000) . This technique is now routinely used to 
study the dark matter distribution of galaxy clusters and 
has recently been detected in the field (Wittman et al 2000; 
van Waerbeke et al 2000; Bacon, Refregier & Elhs 2000; 
Kaiser et al 2000; MaoU et al 2001; Rhodes, Refregier & 
Groth 2001; van Waerbeke et al 2001). AU studies of weak 
lensing have been performed in the optical and IR bands, 
where the images are directly obtained in real space. 
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Fig. 4. — The covariance matrix for the shapelet coefficients of 
nine sources in the FIRST grid pointing. The sources are separated 
by the horizontal and vertical lines. The diagonal entries correspond 
to the variance (i.e., error) of each shapelet coefficient. The block- 
diagonal boxes are the covariance matrix of the coefficients for each 
of the nine sources. The off-diagonal boxes quantify the cross-talk 
between sources. 
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Fig. 5. — The error matrix for one of the FIRST sources plotted in 
the (nl, n2) plane. The shapelet reconstructed image of this source 
is shown on the bottom-left panel of Fig. 3. These variances can also 
be used to Weiner filter the coefficients for image reconstruction, as 
described in 53.4. 



There are a number of reasons to try to extend these 
studies to interferometric images in the radio band. 
Firstly, the brightest radio sources are at high redshift, 
thereby increasing the strength of the lensing signal. Sec- 
ondly, radio interferometers have a well-known and deter- 
ministic convolution beam, and thus do not suffer from 
the irreproducible effects of atmospheric seeing. Thirdly, 
existing surveys such as the FIRST radio survey (Becker 
et al. 1995; White el al. 1997) provide a sparsely sampled 
but very wide-area survey, which offers the unique oppor- 
tunity to measure a weak lensing signal on large angular 



scales (Kamionkowski et al. 1998; Refregier et al. 1998; 
see also Schneider (1999) for the case of SKA). Finally, 
surveys at higher frequencies or in more extended antenna 
configurations could potentially yield very high angular 
resolution and are not limited by the irreducible effects of 
the seeing disk in ground-based optical surveys. 

Because the distortions induced by lensing are only on 
the order of 1%, the shapes of background objects must 
be measured with high precision. In addition, systematic 
effects such as the convolution beam and instrumental dis- 
tortions must be tightly controlled. For this purpose, a 
number of shear measurement methods have been devel- 
oped. The original method of Kaiser, Squires & Broad- 
hurst (1995) was found to be acceptable for current clus- 
ter and large-scale structure surveys (Bacon et al. 2000b; 
Erben et al. 2000), but are not sufficiently reliable for fu- 
ture high-precision surveys. Consequently, several other 
methods have been proposed (Kuijken 1999; Kaiser 2000; 
Rhodes, Refregier & Groth 2000, Berstein & Jarvis 2001). 

Recently, Refregier & Bacon (2001, Paper II) developed 
a new method based on shapelets and demonstrated its 
simplicity and accuracy for ground-based surveys. It is 
thus straightforward to apply this method to interfero- 
metric measurements. Indeed, the shapelet coefficients /„ 
which we derive from the fit on the uv plane (after co- 
adding if required) are already deconvolved from the effec- 
tive beam and can thus be directly used to estimate the 
shear. This can be done using the estimators for the shear 
components 71 and 72 which are given by (see Paper II) 

^' = V s. if V ^^^^ 

where the sum is over even (odd) shapelet coefficients for 
7i (72) E^nd the brackets denote an average over an (un- 
lensed) object ensemble. The matrix Smm is the shear 
matrix, and can be expressed as simple combinations of 
ladder operators in the QHO formalism. These estima- 
tors for individual shapelet components are then optimally 
weighted and combined to provide a minimum-variance 
estimator for the shear. This permits us to achieve the 
highest possible sensitivity (while remaining linear in the 
surface brightness) by using all the available shape infor- 
mation of the lensed sources. 

In Kamionkowski et al. (1998) and Refregier et al. 
(1998), it has been shown that the FIRST radio survey 
is a unique database for measuring weak lensing by large- 
scale structure on large angular scales. In a future paper, 
we will apply the method described here to this survey, 
search for the lensing signal, and, from its amplitude, de- 
rive constraints on cosmological parameters. 

6. CONCLUSIONS 

We have presented a new method for image recon- 
struction from interferometers. Our method is based on 
shapelet decomposition and is simple and robust. It con- 
sists of a linear fit of the shapelet coefficients directly in 
the uv plane, and thus permits a full correction of sys- 
tematic shape distortions caused by bandwidth smearing, 
time-averaging and non-coplanarity Because the fit is lin- 
ear in the shapelet coefficients it can be implemented as 
simple matrix multiplications. It provides the full covari- 
ance matrix of the shapelet coefficients which can then 
be used to estimate errors and cross-talk in the recovered 



shapes of sources. We have shown how source shapes from 
different pointings can be easily co-added using a weighted 
sum of the recovered shapelet coefhcients. We have also 
described how the shapelet parameters could be combined 
to derive optimal image reconstruction, photometry, as- 
trometry and pointing co-addition. 

Our method can be efhciently implemented on paral- 
lel computers. We find that a fit to all the sources in 
a FIRST grid pointing takes about 1 minute on an Ori- 
gin 2000 supercomputer with 10 processors (10 minutes in 
scalar mode). Because we are fitting all sources simulta- 
neously, 0.5 to 1.5 GB of memory is required. 

To test our methods, we considered the observing condi- 
tions of the FIRST radio survey (Becker et al. 1995; White 
et al. 1997) whose snapshot mode yields a sparse sampling 
in uv space. Using numerical simulations tuned to repro- 
duce the conditions of FIRST, we find that the sources 
are well-reconstructed with our method. We have also ap- 
plied our method to a FIRST snapshot pointing and found 
that the shapes are well-recovered. The reconstruction of 
our method compares well with the CLEAN reconstruc- 
tion, without suffering the potential biases inherent in the 
latter method. Our method is thus well-suited for ap- 
plications requiring quantitative and high-precision shape 
measurements . 

In particular, our method is ideal for the measurement 
of the small distortions induced by gravitational lensing in 
the shape of background sources by intervening structures. 
Such a measurement from CLEANed images may well not 
be practical since the systematic distortions induced by 
that method are very difficult to control. (One could per- 
haps imagine running numerical simulations to calibrate 
the shear estimator, but this would be both computation- 
ally expensive and rather uncertain). We have shown how 
our results can be combined with the shear measurement 



method described in Refregier & Bacon (2001) to derive a 
measurement of weak lensing with interferometers. This is 
facilitated both by the fact that our recovered shapelet co- 
efficients are already deconvolved from the effective (dirty) 
beam, and as a consequence of the remarkable properties 
of shapelets under shears. 

Our method therefore opens the possibility of high- 
precision measurements of weak lensing with interferom- 
eters. While to date all weak-lensing studies have been 
carried using optical data (and therefore in real space) , an 
interferometric measurement of weak lensing in the radio 
band is very attractive (Kamionkowski et al. 1998; Re- 
fregier et al. 1998; Schneider 1999). Indeed, the lensing 
signal is expected to be larger because radio sources have 
a higher mean redshift. In addition, such a measurement 
would not suffer from the irreproducible effects of atmo- 
spheric seeing. Instead, the effective (dirty) beam is fully 
known for interferometers and the noise properties of the 
antennas are well-understood. As a result, the impact of 
systematic effects, the crucial limitation in the search for 
weak lensing, are expected to be lower with radio inter- 
ferometers. In a future paper, we will describe our mea- 
surement of weak lensing by large-scale structure with the 
FIRST survey using the present method. 
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